{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "736e5076",
   "metadata": {},
   "source": [
    "# Chapter 3 - Optimal Flows (Julia Code)\n",
    "\n",
    "## Bellman’s Method\n",
    "\n",
    "Here we demonstrate solving a shortest path problem using Bellman's method.\n",
    "\n",
    "Our first step is to set up the cost function, which we store as an array\n",
    "called `c`. Note that we set `c[i, j] = Inf` when no edge exists from `i` to\n",
    "`j`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c36f47c4",
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "c = fill(Inf, (7, 7))\n",
    "c[1, 2], c[1, 3], c[1, 4] = 1, 5, 3\n",
    "c[2, 4], c[2, 5] = 9, 6\n",
    "c[3, 6] = 2\n",
    "c[4, 6] = 4\n",
    "c[5, 7] = 4\n",
    "c[6, 7] = 1\n",
    "c[7, 7] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b8fa3ad",
   "metadata": {},
   "source": [
    "Next we define the Bellman operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2bc51e6",
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "function T(q)\n",
    "    Tq = similar(q)\n",
    "    n = length(q)\n",
    "    for x in 1:n\n",
    "        Tq[x] = minimum(c[x, :] + q[:])\n",
    "    end\n",
    "    return Tq\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9857bf8a",
   "metadata": {},
   "source": [
    "Now we arbitrarily set $q \\equiv 0$, generate the sequence of iterates $T_q$,\n",
    "$T^2_q$, $T^3_q$ and plot them. By $T^3_q$ has already converged on $q^∗$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "757cadc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyPlot\n",
    "export_figures = false\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "\n",
    "n = 7\n",
    "q = zeros(n)\n",
    "ax.plot(1:n, q)\n",
    "ax.set_xlabel(\"cost-to-go\")\n",
    "ax.set_ylabel(\"nodes\")\n",
    "\n",
    "for i in 1:3\n",
    "    new_q = T(q)\n",
    "    ax.plot(1:n, new_q, \"-o\", alpha=0.7, label=\"iterate $i\")\n",
    "    q = new_q\n",
    "end\n",
    "\n",
    "ax.legend()\n",
    "if export_figures == true\n",
    "    plt.savefig(\"figures/shortest_path_iter_1.pdf\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d31abe43",
   "metadata": {},
   "source": [
    "## Linear programming\n",
    "\n",
    "When solving linear programs, one option is to use a domain specific modeling\n",
    "language to set out the objective and constraints in the optimization problem.\n",
    "Here we demonstrate the Julia package `JuMP`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67361f35",
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuMP\n",
    "using GLPK"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bd4aef4",
   "metadata": {},
   "source": [
    "We create our model object and select our solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "856ef0ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Model()\n",
    "set_optimizer(m, GLPK.Optimizer)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe817ce7",
   "metadata": {},
   "source": [
    "Now we add variables, constraints and an objective to our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07228cb0",
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "@variable(m, q1 >= 0)\n",
    "@variable(m, q2 >= 0)\n",
    "@constraint(m, 2q1 + 5q2 <= 30)\n",
    "@constraint(m, 4q1 + 2q2 <= 20)\n",
    "@objective(m, Max, 3q1 + 4q2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "842781b0",
   "metadata": {},
   "source": [
    "Finally we solve our linear program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83c0d928",
   "metadata": {},
   "outputs": [],
   "source": [
    "optimize!(m)\n",
    "\n",
    "println(value.(q1)) \n",
    "println(value.(q2))"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "Julia",
   "language": "Julia",
   "name": "julia-1.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
